## Compositional time series analysis of state budgets
## using seemingly unrelated regression of additive logratios
## with optional lags and state fixed effects
## Christopher Adolph, Christian Breunig, and Chris Koski
## 3 March 2018
##                                        
## Model 4:  Control for TELs
##
## This program does the following:
##
## 1. Loads compositional data and covariates
## 2. Runs a SUR model of those data, with optional lags and fixed effects
##    (Specifically, the SUR models the additive logratios of the composition.
##    All equations have the same covariates, and only lags of the current
##    logratio are included in a given equation.)
## 3. Simulates expected compositions under different values of the
##    covariates
##    (Specifically, the composition is assumed to rest initially at its
##    sample mean.   All covariates are initially at their sample means.
##    We then change a single covariate or combination of covariates,
##    and iterate the model forwards in time under the new
##    "counterfactual scenario".)
## 4. Plots relative ratios of compositions by scenario of covariate values
#
## User input is required first to set up the model (the first user options
## block) and to set up the counterfactuals (the second user options block).
##          

# Clear memory
rm(list = ls())

#######################################################################
# User settings
# (Create a new file for each composition + specification)
dsname <- "dataABK.csv"    # Must be comma-seperated variable file
specname <- "M4"

lhs <- c("K12Ed",
         "Welfare",
         "Medical",         
         "NaturalResources",
         "HigherEd",
         "Highways",
         "PublicSafety",
         "Miscellaneous"
         )

lhsNice <- c("K-12\nEducation",
             "Medicaid\n& Welfare",
             "Public Health\n& Hospitals",
             "Natural\nResources",
             "Higher\nEducation",
             "Highways",
             "Police\n& Prisons",
             "Other\nSpending"
             )

lhsWide <- c("K-12",
             "Medicaid & Welfare",
             "Public Health & Hospitals",
             "Natural Resources",
             "Higher Ed",
             "Highways",
             "Police & Prisons",
             "Other Spending"
             )

## Covariates
rhs <- c("ACIR", 
         "govadd",
         "DUG", "RUG",
         "rpincpc",
         "popdensity",
         "Unemployment",
         "South","Northeast","West",
         "trend",
         "pctunder19",
         "pctover65",
         "TEL"
         )
                              # Covariates
                              # Note: interactions must be constructed in
                              # advance (i.e., exist as variables in the
                              # datafile)
numlags <- 1                  # Number of lags in model
csfe <- 0                     # Include state fixed effects? (0/1)
csvar <- "State"              # Name of state variable
csvarid <- "FIPS"             # Name of state id number variable
tsvar <- "Year"               # Name of time variable
sims <- 100000                # Num of simulations for post-estimation
ci <- c(0.95, 0.90)           # Confidence intervals of first diffs
simperiods <- 4               # Periods predicted ahead in simulations
plotperiod <- simperiods      # Period plotted in ropeladders
xlagnum <- 1                  # Lag for covariates
# User must also supply counterfactuals below
#######################################################################

# Load libraries (make sure these are available to R)
library(tile)
library(simcf)
# Download tile and simcf from:
#  http://faculty.washington.edu/cadolph/software
library(compositions)
library(MASS)
library(systemfit)
source("helper.r")     # Keep this file in working directory

# Load data
data <- read.csv(dsname, na.strings = "NA", header=TRUE)

## Delete early years
data <- data[data$Year>=1983,]

## Delete NAs in 2010
data <- data[data$Year<2010,]

## Remove Alaska and Hawaii
data <- data[!(data$State=="AK")&!(data$State=="HI")&!(data$State=="NE"),]

attach(data)
c <- getvar(csvarid)
t <- getvar(tsvar)

# Set endogenous variables
components <- getvar(lhs)       # need to fix var names
components <- clo(components)   # close the composition
alrcomp <- alr(components)      # convert to centered logratios
alrcomp.names <- names(alrcomp)
for (i in 1:length(alrcomp.names))
  alrcomp.names[i] <- paste(alrcomp.names[i],"LR",sep="")
names(alrcomp) <- alrcomp.names
datarun <- alrcomp
ncomp <- length(lhs)
compavg <- apply(na.omit(components),2,"mean")

# Make lags of response variables
lagcomp <- lagcomp.names <- NULL
for (i in 1:ncol(alrcomp)) {
  for (j in 1:numlags) {
    curlag <- paste("lag",j,sep="")
    lagcomp <- cbind(lagcomp,lagpanel(alrcomp[,i],c,t,j))
    lagcomp.names <- c(lagcomp.names,paste(names(alrcomp)[i],curlag,sep=""))
  }
}
colnames(lagcomp) <- lagcomp.names
datarun <- cbind(datarun,as.data.frame(lagcomp))

# Add covariates and constant
xcovariates <- getvar(rhs)
lagcovariates <- as.data.frame(lagpanel(xcovariates,c,t,xlagnum))
names(lagcovariates) <- names(xcovariates)
xcovariates <- lagcovariates
datarun <- cbind(datarun,xcovariates)
constant <- as.data.frame(matrix(data=1,nrow=nrow(datarun),ncol=1))
names(constant) <- "constant"
datarun <- cbind(datarun,constant)
numcovariates <- length(rhs)

# Add fixed effects
if (csfe) {
  stateid <- getvar(csvar)
  cslist <- unique(stateid)
  ncs <- nrow(cslist)
  dummies <- matrix(nrow=nrow(datarun),ncol=nrow(cslist))
  for (i in 1:nrow(datarun))
    dummies[i,] <- as.numeric(as.vector(stateid[i,1])==as.vector(cslist))
  colnames(dummies) <- as.vector(t(cslist))
  datarun <- cbind(datarun,as.data.frame(dummies))
} else {
  cslist <- "constant"
}

# Construct system of equations
eq <- NULL
for (i in 1:ncomp-1) {
  eq[i] <- paste(names(datarun)[i],"~","- 1")
  for (j in 1:numlags) {
    curlagnames <- lagcomp.names[(i-1)*numlags+1:(i*numlags)]
    eq[i] <- paste(eq[i],"+",curlagnames[j])
  }
  for (j in 1:numcovariates) {
    eq[i] <- paste(eq[i],"+",rhs[j])
  }
  if (csfe) {
    for (k in 1:ncs) {
      eq[i] <- paste(eq[i],"+",cslist[k,1])
    }
  } else {
    eq[i] <- paste(eq[i],"+","constant")
  }
}
eqns <- list(comp1="y~x")
eqnnames <- NULL
ieq <- 1
while(ieq<ncomp) {
  ieqn <- paste("comp",ieq,sep="")
  eqnnames <- c(eqnnames,ieqn)
  cureq <- eq[ieq]
  eqns[[ieqn]] <- cureq
  eqns[[ieqn]] <- as.formula(eqns[[ieqn]])
  ieq <- ieq+1
}


# Listwise deletion of missings
datarun <- na.omit(datarun)

# Run SUR
sur.res1 <- systemfit("SUR",
                      formula=eqns,
                      data=datarun,
                      control = systemfit.control(maxiter=100)
                      )
print(summary(sur.res1))

# Extract b and vcov
b <- sur.res1$coefficients
vc <- sur.res1$coefCov
write.csv(printnicetable(b,vc,7),paste(specname,".csv",sep=""))

# Draw from predictive distribution of parameters
simbetas <- mvrnorm(n = sims, mu=b, Sigma=vc)


#############################################################################
## Set counterfactuals:  User supplied
nscen <- 17              # Specify the number of counterfactual scenarios run
curscen <- 0

# Set up a template for all scenarios (sets all covariates at means)
xscen <- setx.make(nscen,xcovariates,as.vector(t(cslist)))
# set up scenario 1
curscen <- curscen + 1
xscen <- setx.name(xscen,"Rep -> Unified Dem Control",scen=curscen)
xscen <- setx.change(xscen,"DUG",1,0,scen=curscen)
xscen <- setx.change(xscen,"RUG",0,1,scen=curscen)

# set up scenario 2
curscen <- curscen + 1
xscen <- setx.name(xscen,"Div -> Unified Dem Control",scen=curscen)
xscen <- setx.change(xscen,"DUG",1,0,scen=curscen)
xscen <- setx.change(xscen,"RUG",0,0,scen=curscen)

# set up scenario 3
curscen <- curscen + 1
xscen <- setx.name(xscen,"Div -> Unified Rep Control",scen=curscen)
xscen <- setx.change(xscen,"RUG",1,0,scen=curscen)
xscen <- setx.change(xscen,"DUG",0,0,scen=curscen)

# set up scenario 4
curscen <- curscen + 1
xscen <- setx.name(xscen,"Dem -> Unified Rep Control",scen=curscen)
xscen <- setx.change(xscen,"RUG",1,0,scen=curscen)
xscen <- setx.change(xscen,"DUG",0,1,scen=curscen)

# set up scenario 5
curscen <- curscen + 1
xscen <- setx.name(xscen,"Powerful Governor",scen=curscen)
xscen <- setx.change(xscen,"govadd",mean(govadd) + sd(govadd), mean(govadd),scen=curscen)

# set up scenario 6
curscen <- curscen + 1
xscen <- setx.name(xscen,"Stringent Budget",scen=curscen)
xscen <- setx.change(xscen,"ACIR", max(ACIR), mean(ACIR),scen=curscen)

# set up scenario 7
curscen <- curscen + 1
xscen <- setx.name(xscen,"Tax or Expenditure\nLimit",scen=curscen)
xscen <- setx.change(xscen,"TEL", 1, 0, scen=curscen)

# Set up scenario 8
curscen <- curscen + 1
xscen <- setx.name(xscen,"Unemployment\nRate",scen=curscen)
xscen <- setx.change(xscen,"Unemployment",mean(Unemployment)+ sd(Unemployment),mean(Unemployment),scen=curscen)

# set up scenario 9
curscen <- curscen + 1
xscen <- setx.name(xscen,"Real Income\nPer Capita",scen=curscen)
xscen <- setx.change(xscen,"rpincpc",mean(rpincpc) + sd(rpincpc),mean(rpincpc),scen=curscen)

# set up scenario 10
curscen <- curscen + 1
xscen <- setx.name(xscen,"Population\nDensity",scen=curscen)
xscen <- setx.change(xscen,"popdensity",mean(popdensity) + sd(popdensity), mean(popdensity),scen=curscen)

# Set up scenario 11
curscen <- curscen + 1
xscen <- setx.name(xscen,"Young\n(<20 Years)",scen=curscen)
xscen <- setx.change(xscen,"pctunder19",mean(pctunder19)+sd(pctunder19),
                     mean(pctunder19),scen=curscen)
rpcfVal <- rpcf(c(mean(pctunder19), mean(pctover65), 1- mean(pctunder19) - mean(pctover65)),
                1, delta=sd(pctunder19))
xscen <- setx.change(xscen,"pctover65",rpcfVal[2],mean(pctover65),scen=curscen)

# Set up scenario 12
curscen <- curscen + 1
xscen <- setx.name(xscen,"Adults\n(20 to 64 Years)",scen=curscen)
midageMN <- mean(1- pctunder19 -pctover65)
midageSD <- sd(1- pctunder19 -pctover65) 
rpcfVal <- rpcf(c(mean(pctover65), mean(pctunder19), midageMN),
                3, delta=midageSD)
xscen <- setx.change(xscen,"pctover65",rpcfVal[1],mean(pctover65),scen=curscen)
xscen <- setx.change(xscen,"pctunder19",rpcfVal[2],mean(pctunder19),scen=curscen)

# Set up scenario 13
curscen <- curscen + 1
xscen <- setx.name(xscen,"Elderly\n(>=65 Years)",scen=curscen)
xscen <- setx.change(xscen,"pctover65",mean(pctover65)+sd(pctover65),mean(pctover65),scen=curscen)
rpcfVal <- rpcf(c(mean(pctover65), mean(pctunder19), 1- mean(pctunder19) - mean(pctover65)),
                1, delta=sd(pctover65))
xscen <- setx.change(xscen,"pctunder19",rpcfVal[2],mean(pctunder19),scen=curscen)

# Set up scenario 14
curscen <- curscen + 1
xscen <- setx.name(xscen,"Northeast",scen=curscen)
denom <- 1
xscen <- setx.change(xscen,"Northeast",1,0,scen=curscen)
xscen <- setx.change(xscen,"West",0,mean(na.omit(West))/denom,scen=curscen)
xscen <- setx.change(xscen,"South",0,mean(na.omit(South))/denom,scen=curscen)

# Set up scenario 15
curscen <- curscen + 1
xscen <- setx.name(xscen,"Midwest",scen=curscen)
denom <- 1
xscen <- setx.change(xscen,"Northeast",0,mean(na.omit(Northeast))/denom,scen=curscen)
xscen <- setx.change(xscen,"South",0,mean(na.omit(South))/denom,scen=curscen)
xscen <- setx.change(xscen,"West",0,mean(na.omit(West))/denom,scen=curscen)

# Set up scenario 16
curscen <- curscen + 1
xscen <- setx.name(xscen,"South",scen=curscen)
denom <- 1
xscen <- setx.change(xscen,"South",1,0,scen=curscen)
xscen <- setx.change(xscen,"Northeast",0,mean(na.omit(Northeast))/denom,scen=curscen)
xscen <- setx.change(xscen,"West",0,mean(na.omit(West))/denom,scen=curscen)

# Set up scenario 17
curscen <- curscen + 1
xscen <- setx.name(xscen,"West",scen=curscen)
denom <- 1
xscen <- setx.change(xscen,"West",1,0,scen=curscen)
xscen <- setx.change(xscen,"South",0,mean(na.omit(South))/denom,scen=curscen)
xscen <- setx.change(xscen,"Northeast",0,mean(na.omit(Northeast))/denom,scen=curscen)



#############################################################################
# Loop over counterfactuals
simestall <- simlowall <- simupall <- simseall <- NULL
simfdestall <- simfdlowall <- simfdupall <- simfdseall <- NULL
simfrestall <- simfrlowall <- simfrupall <- simfrlowall0 <- simfrupall0 <- simfrseall <- NULL
scendescriptall <- NULL
for (iscen in 1:nrow(xscen$xhyp)) {

  # Simulate dynamics from estimated model
  nparam <- NULL
  for (i in 1:ncomp-1)
    nparam[i] <- length(sur.res1$eq[[as.character(i)]]$coefficients)
  simest <- simlow <- simup <- simse <- matrix(nrow=simperiods,ncol=ncomp)
  simfdest <- simfdlow <- simfdup <- simfdse <- matrix(nrow=simperiods,ncol=ncomp)
  simfrest <- simfrlow <- simfrup <- simfrlow0 <- simfrup0 <- simfrse <- matrix(nrow=simperiods,ncol=ncomp)
  simy <- simy0 <- list()
  for (i in 1:simperiods) {
    simy[[i]] <- matrix(nrow=sims,ncol=(ncomp-1))
    simy0[[i]] <- matrix(nrow=sims,ncol=(ncomp-1))
  }
  basecol <- 0
  ieq <- 1
  while (ieq<ncomp) {

    xhyp <- c(rep(mean(alrcomp[,ieq]),numlags), xscen$xhyp[iscen,])
    xhyp0 <- c(rep(mean(alrcomp[,ieq]),numlags), xscen$xhyp0[iscen,])

    ieqn <- paste("comp",ieq,sep="")
    simbetascur <- simbetas[,(basecol+1):(basecol+length(sur.res1$eq[[ieq]]$coefficients))]
    simcomp <- postestd(xhyp=xhyp,
                        xhyp0=xhyp0,
                        t=simperiods,
                        sims=sims,
                        ci=ci[1],
                        lagcols=seq(from=1, to=numlags),
                        simbetas=simbetascur,
                        dofirstdiff=FALSE
                        )

    for (iperiod in 1:simperiods) {
      simy[[iperiod]][,ieq] <- simcomp$simy[,iperiod]
      simy0[[iperiod]][,ieq] <- simcomp$simy0[,iperiod]
    }
    basecol <- basecol + length(sur.res1$eq[[ieq]]$coefficients)
    ieq <- ieq+1
  }

  # Transform results to composition space & construct ci's
  for (iperiod in 1:simperiods) {
    simy[[iperiod]] <- alrInv(simy[[iperiod]])
    simy0[[iperiod]] <- alrInv(simy0[[iperiod]])
    for (icomp in 1:ncomp) {
      simysort <- sort(simy[[iperiod]][,icomp])
      simest[iperiod,icomp] <- mean(simysort)
      simse[iperiod,icomp] <- sd(simysort)
      simup[iperiod,icomp] <- simysort[trunc(sims*((1+ci[1])/2))]
      simlow[iperiod,icomp] <- simysort[trunc(sims*((1-ci[1])/2))]

      simfdsort <- sort(simy[[iperiod]][,icomp] - simy0[[iperiod]][,icomp])
      simfdest[iperiod,icomp] <- mean(simfdsort)
      simfdse[iperiod,icomp] <- sd(simfdsort)
      simfdup[iperiod,icomp] <- simfdsort[trunc(sims*((1+ci[1])/2))]
      simfdlow[iperiod,icomp] <- simfdsort[trunc(sims*((1-ci[1])/2))]

      simfrsort <- sort(simy[[iperiod]][,icomp] / simy0[[iperiod]][,icomp])
      simfrest[iperiod,icomp] <- mean(simfrsort)
      simfrse[iperiod,icomp] <- sd(simfrsort)
      simfrup[iperiod,icomp] <- simfrsort[trunc(sims*((1+ci[1])/2))]
      simfrlow[iperiod,icomp] <- simfrsort[trunc(sims*((1-ci[1])/2))]
      simfrup0[iperiod,icomp] <- simfrsort[trunc(sims*((1+ci[2])/2))]
      simfrlow0[iperiod,icomp] <- simfrsort[trunc(sims*((1-ci[2])/2))]
    }
  }

  simestall <- rbind(simestall,simest[plotperiod,])
  simlowall <- rbind(simlowall,simlow[plotperiod,])
  simupall <- rbind(simupall,simup[plotperiod,])
  scendescriptall <- c(scendescriptall,xscen$name[iscen])

  simfdestall <- rbind(simfdestall,simfdest[plotperiod,])
  simfdlowall <- rbind(simfdlowall,simfdlow[plotperiod,])
  simfdupall <- rbind(simfdupall,simfdup[plotperiod,])
    
  simfrestall <- rbind(simfrestall,simfrest[plotperiod,])
  simfrlowall <- rbind(simfrlowall,simfrlow[plotperiod,])
  simfrupall <- rbind(simfrupall,simfrup[plotperiod,])
  simfrlowall0 <- rbind(simfrlowall0,simfrlow0[plotperiod,])
  simfrupall0 <- rbind(simfrupall0,simfrup0[plotperiod,])
}




### Make Figure 2 

peALL <- simfrestall[1:4,]
lwALL <- simfrlowall[1:4,]
upALL <- simfrupall[1:4,]
lwALL0 <- simfrlowall0[1:4,]
upALL0 <- simfrupall0[1:4,]

colnames(peALL) <- lhsNice

row.names(peALL) <- c("REP -> DEM",
                      "DIV -> DEM",
                      "DIV -> REP",
                      "DEM -> REP"                      
                      )

collectedTraces <- collectedBlanks <- vector("list", nrow(peALL))
for (i in 1:nrow(peALL)) {

    labelscur <- row.names(peALL)
        
    pch <- rep(1,ncol(peALL))
    for (j in 1:ncol(peALL))
        if (((lwALL[i,j]-1)*(upALL[i,j]-1))>0) {
            pch[j] <- 16
        } 
    
    collectedTraces[[i]] <-  ropeladder(x=peALL[i,],
                                        lower=lwALL[i,],
                                        upper=upALL[i,],
                                        labels=colnames(peALL),
                                        sublabels=labelscur[i],
                                        sublabelsX=0.1,
                                        entryheight=0.175,
                                        subentryheight=1.0,
                                        size=0.6,
                                        pch=pch,
                                        col="black",
                                        spaceAbove=0.05,
                                        plot=1
                                        )
}

text1 <- textTile(x=0.685,
                  y=0.965,
                  labels="Change in budget",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)

text2 <- textTile(x=0.795,
                  y=0.965,
                  labels="if control shifts...",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)       

vertmark <- linesTile(x=c(1,1), y=c(0,1), plot=1)


repPrioritiesLab <- textTile(labels="REPUBLICAN PRIORITIES",
                             x=0.62,
                             y=0.2385,
                             rot=90,
                             cex=0.75,
                             col="black",
                             clip="off",
                             alpha=1,
                             plot=1)

repPrioritiesBox <- polygonTile(x=c(0.61, 0.61, 0.63, 0.63),
                                y=c(0.0195, 0.4575, 0.4575, 0.0195),
                                fill="gray75",
                                col=NA,
                                clip="off",
                                plot=1)


demPrioritiesLab <- textTile(labels="DEM PRIORITIES",
                             x=0.62,
                             y=0.83,
                             rot=90,
                             cex=0.75,
                             col="white",
                             clip="off",
                             alpha=1,
                             plot=1)

demPrioritiesBox <- polygonTile(x=c(0.61, 0.61, 0.63, 0.63),
                                y=c(0.73, 0.93, 0.93, 0.73),
                                fill="gray35",
                                col=NA,
                                clip="off",
                                plot=1)

limits <- c(0.75,1.175)
at <- c(0.85,0.9,0.95,1,1.05,1.1,1.15)
labels <- c("-15%","-10%","-5%", "0%", "+5%", "+10%","+15%")

# Create Plot 2 and save to pdf -- first pass
file <- paste0("partisanBudgetEffects", specname)
tc <- tile(collectedTraces,
     vertmark, text1, text2,
     repPrioritiesLab, repPrioritiesBox,
     demPrioritiesLab, demPrioritiesBox,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file, width=6.5),
     special=list(fontsizeRopeladder=5.5)
     )
           
ycoords <- NULL
for (i in 1:4) ycoords <- rbind(ycoords, tc$traces[[i]]$y)

peALL0 <- peALL
pch <- peALL0*0 + 16
col <- matrix("white", nrow=nrow(peALL), ncol=ncol(peALL))
for (j in 1:ncol(peALL)) {
    for (i in 1:nrow(peALL)) {
        sig95 <- ((lwALL[i,j]-1)*(upALL[i,j]-1))>0
        sig90 <- ((lwALL0[i,j]-1)*(upALL0[i,j]-1))>0
        
        if (sig95&sig90) {
            peALL0[i,j] <- NA
        } else {
            if (sig90) {
                col[i,j] <- "gray65"
            }
        }
    }
}

nsTrace <- pointsTile(x=as.numeric(peALL0),
                      y=as.numeric(ycoords),
                      col=as.character(col),
                      pch=pch,
                      alpha=1,
                      layer=1,
                      size=0.51,
                      plot=1)

# Create Plot 2 and save to pdf -- second pass
file <- paste0("partisanBudgetEffects", specname)
tile(collectedTraces, nsTrace,
     vertmark, text1, text2,
     repPrioritiesLab, repPrioritiesBox,
     demPrioritiesLab, demPrioritiesBox,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file, width=6.5),
     special=list(fontsizeRopeladder=5.5)
     )



## Make Plot 3

peALL <- simfrestall[5:7,,drop=FALSE]
lwALL <- simfrlowall[5:7,,drop=FALSE]
upALL <- simfrupall[5:7,,drop=FALSE]
lwALL0 <- simfrlowall0[5:7,,drop=FALSE]
upALL0 <- simfrupall0[5:7,,drop=FALSE]

colnames(peALL) <- lhsWide

row.names(peALL) <- c("Power of\nGovernor\n(Beyle)",
                      "Budget\nStringency",
                      "Tax or\nExpenditure\nLimit"
                      )

collectedTraces <- collectedBlanks <- vector("list", nrow(peALL))
for (i in 1:ncol(peALL)) {

    labelscur <- colnames(peALL)
    
    pch <- rep(1,nrow(peALL))
    for (j in 1:nrow(peALL))
        if (((lwALL[j,i]-1)*(upALL[j,i]-1))>0) {
            pch[j] <- 16
        }
    
    collectedTraces[[i]] <-  ropeladder(x=peALL[,i],
                                        lower=lwALL[,i],
                                        upper=upALL[,i],
                                        labels=row.names(peALL),
                                        sublabels=labelscur[i],
                                        sublabelsX=0.135,
                                        entryheight=0.400,
                                        subentryheight=1.20,
                                        size=0.6,
                                        pch=pch,
                                        spaceAbove=0.025,
                                        plot=1
                                        )
}

text1 <- textTile(x=0.694,
                  y=0.97,
                  labels="+1 sd in...",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)

text2 <- textTile(x=0.8075,
                  y=0.97,
                  labels="shifts these budgets",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)
                  
vertmark <- linesTile(x=c(1,1), y=c(0,1), plot=1)

limits <- c(0.75,1.175)
at <- c(0.85,0.9,0.95,1,1.05,1.1,1.15)
labels <- c("-15%","-10%","-5%", "0%", "+5%", "+10%","+15%")

# Create plot and save to pdf -- first pass
file <- paste0("instBudgetEffects", specname)
tc <- tile(collectedTraces,
     vertmark, text1, text2,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file,width=6.5),
     special=list(fontsizeRopeladder=5.5)
     )

ycoords <- NULL
for (i in 1:length(lhs)) ycoords <- rbind(ycoords, tc$traces[[i]]$y)

peALL0 <- peALL
pch <- peALL0*0 + 16
col <- matrix("white", nrow=nrow(peALL), ncol=ncol(peALL))
for (i in 1:ncol(peALL)) {
    for (j in 1:nrow(peALL)) {
        sig95 <- ((lwALL[j,i]-1)*(upALL[j,i]-1))>0
        sig90 <- ((lwALL0[j,i]-1)*(upALL0[j,i]-1))>0
        
        if (sig95&sig90) {
            peALL0[j,i] <- NA
        } else {
            if (sig90) {
                col[j,i] <- "gray65"
            }
        }
    }
}

nsTrace <- pointsTile(x=as.numeric(t(peALL0)),
                      y=as.numeric(ycoords),
                      col=as.character(t(col)),
                      pch=pch,
                      alpha=1,
                      layer=1,
                      size=0.51,
                      plot=1)

# Create plot and save to pdf -- second pass
file <- paste0("instBudgetEffects", specname)
tile(collectedTraces, nsTrace,
     vertmark, text1, text2,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file, width=6.5),
     special=list(fontsizeRopeladder=5.5)
     )


## Make Plot 4

peALL <- simfrestall[8:9,]
lwALL <- simfrlowall[8:9,]
upALL <- simfrupall[8:9,]
lwALL0 <- simfrlowall0[8:9,]
upALL0 <- simfrupall0[8:9,]

colnames(peALL) <- lhsWide

row.names(peALL) <- c("Unemployment\nRate",
                      "Real Income\nPer Capita"
                      )

collectedTraces <- collectedBlanks <- vector("list", nrow(peALL))
for (i in 1:ncol(peALL)) {

    labelscur <- colnames(peALL)
    
    pch <- rep(1,nrow(peALL))
    for (j in 1:nrow(peALL))
        if (((lwALL[j,i]-1)*(upALL[j,i]-1))>0) {
            pch[j] <- 16
        }
    
    collectedTraces[[i]] <-  ropeladder(x=peALL[,i],
                                        lower=lwALL[,i],
                                        upper=upALL[,i],
                                        labels=row.names(peALL),
                                        sublabels=labelscur[i],
                                        sublabelsX=0.135,
                                        entryheight=0.425,
                                        subentryheight=1.35,
                                        size=0.6,
                                        pch=pch,
                                        spaceAbove=0.05,
                                        plot=1
                                        )
}

text1 <- textTile(x=0.675,
                  y=0.955,
                  labels="+1 sd in...",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)

text2 <- textTile(x=0.8075,
                  y=0.955,
                  labels="shifts these budgets",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)
                  
vertmark <- linesTile(x=c(1,1), y=c(0,1), plot=1)

limits <- c(0.75,1.175)
at <- c(0.85,0.9,0.95,1,1.05,1.1,1.15)
labels <- c("-15%","-10%","-5%", "0%", "+5%", "+10%","+15%")

# Create plot and save to pdf -- first pass
file <- paste0("econBudgetEffects", specname)
tc <- tile(collectedTraces,
     vertmark, text1, text2,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file, width=6.5),
     special=list(fontsizeRopeladder=5.5)
     )

ycoords <- NULL
for (i in 1:length(lhs)) ycoords <- rbind(ycoords, tc$traces[[i]]$y)

peALL0 <- peALL
pch <- peALL0*0 + 16
col <- matrix("white", nrow=nrow(peALL), ncol=ncol(peALL))
for (i in 1:ncol(peALL)) {
    for (j in 1:nrow(peALL)) {
        sig95 <- ((lwALL[j,i]-1)*(upALL[j,i]-1))>0
        sig90 <- ((lwALL0[j,i]-1)*(upALL0[j,i]-1))>0
        
        if (sig95&sig90) {
            peALL0[j,i] <- NA
        } else {
            if (sig90) {
                col[j,i] <- "gray65"
            }
        }
    }
}

nsTrace <- pointsTile(x=as.numeric(t(peALL0)),
                      y=as.numeric(ycoords),
                      col=as.character(t(col)),
                      pch=pch,
                      alpha=1,
                      layer=1,
                      size=0.51,
                      plot=1)

# Create plot and save to pdf -- second pass
file <- paste0("econBudgetEffects", specname)
tile(collectedTraces, nsTrace,
     vertmark, text1, text2,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file,width=6.5,family="GillSans"),
     special=list(fontsizeRopeladder=5.5)
     )


## Make Plot 4 long version

peALL <- rbind(simfrestall[8:9,], rep(NA, 8))
lwALL <- rbind(simfrlowall[8:9,], rep(NA, 8))
upALL <- rbind(simfrupall[8:9,], rep(NA, 8))
lwALL0 <- rbind(simfrlowall0[8:9,], rep(NA, 8))
upALL0 <- rbind(simfrupall0[8:9,], rep(NA, 8))

colnames(peALL) <- lhsWide

row.names(peALL) <- c("Unemployment\nRate",
                      "Real Income\nPer Capita",
                      "Real Growth in\nTotal Spending"
                      )

collectedTraces <- collectedBlanks <- vector("list", nrow(peALL))
for (i in 1:ncol(peALL)) {

    labelscur <- colnames(peALL)
    
    pch <- rep(1,nrow(peALL))
    for (j in 1:nrow(peALL))
        if (!is.na(lwALL[j,i]))
            if (((lwALL[j,i]-1)*(upALL[j,i]-1))>0) {
                pch[j] <- 16
            }
    
     collectedTraces[[i]] <-  ropeladder(x=peALL[,i],
                                        lower=lwALL[,i],
                                        upper=upALL[,i],
                                        labels=row.names(peALL),
                                        sublabels=labelscur[i],
                                        sublabelsX=0.135,
                                        entryheight=0.400,
                                        subentryheight=1.20,
                                        size=0.6,
                                        pch=pch,
                                        spaceAbove=0.025,
                                        plot=1
                                        )
}

text1 <- textTile(x=0.675,
                  y=0.97,
                  labels="+1 sd in...",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)

text2 <- textTile(x=0.8075,
                  y=0.97,
                  labels="shifts these budgets",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)
                  
vertmark <- linesTile(x=c(1,1), y=c(0,1), plot=1)

limits <- c(0.75,1.175)
at <- c(0.85,0.9,0.95,1,1.05,1.1,1.15)
labels <- c("-15%","-10%","-5%", "0%", "+5%", "+10%","+15%")

# Create plot and save to pdf -- first pass
file <- paste0("econBudgetEffectsLong", specname)
tc <- tile(collectedTraces,
     vertmark, text1, text2,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file, width=6.5),
     special=list(fontsizeRopeladder=5.5)
     )

ycoords <- NULL
for (i in 1:length(lhs)) ycoords <- rbind(ycoords, tc$traces[[i]]$y)

peALL0 <- peALL
pch <- peALL0*0 + 16
col <- matrix("white", nrow=nrow(peALL), ncol=ncol(peALL))
for (i in 1:ncol(peALL)) {
    for (j in 1:nrow(peALL)) {
        if (!is.na(lwALL[j,i])) {
            sig95 <- ((lwALL[j,i]-1)*(upALL[j,i]-1))>0
            sig90 <- ((lwALL0[j,i]-1)*(upALL0[j,i]-1))>0
        } else {
            sig95 <- sig90 <- TRUE
        }
            
        if (sig95&sig90) {
            peALL0[j,i] <- NA
        } else {
            if (sig90) {
                col[j,i] <- "gray65"
            }
        }
    }
}

nsTrace <- pointsTile(x=as.numeric(t(peALL0)),
                      y=as.numeric(ycoords),
                      col=as.character(t(col)),
                      pch=16,
                      alpha=1,
                      layer=1,
                      size=0.51,
                      plot=1)

# Create plot and save to pdf -- second pass
file <- paste0("econBudgetEffectsLong", specname)
tile(collectedTraces, nsTrace,
     vertmark, text1, text2,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file, width=6.5),
     special=list(fontsizeRopeladder=5.5)
     )



## Make Plot 5

peALL <- simfrestall[10:13,]
lwALL <- simfrlowall[10:13,]
upALL <- simfrupall[10:13,]
lwALL0 <- simfrlowall0[10:13,]
upALL0 <- simfrupall0[10:13,]

colnames(peALL) <- lhsWide

row.names(peALL) <- c("Population\nDensity",
                      "Under 19\nYear Olds",
                      "20 to 64\nYear Olds",
                      "65+\nYear Olds"
                      )

collectedTraces <- collectedBlanks <- vector("list", nrow(peALL))
for (i in 1:ncol(peALL)) {

    labelscur <- colnames(peALL)
    
    pch <- rep(1,nrow(peALL))
    for (j in 1:nrow(peALL))
        if (((lwALL[j,i]-1)*(upALL[j,i]-1))>0) {
            pch[j] <- 16
        }
    
    collectedTraces[[i]] <-  ropeladder(x=peALL[,i],
                                        lower=lwALL[,i],
                                        upper=upALL[,i],
                                        labels=row.names(peALL),
                                        sublabels=labelscur[i],
                                        sublabelsX=0.135,
                                        entryheight=0.375,
                                        subentryheight=1.10,
                                        size=0.6,
                                        pch=pch,
                                        spaceAbove=0.025,
                                        plot=1
                                        )
}

text1 <- textTile(x=0.70,
                  y=0.975,
                  labels="+1 sd in...",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)

text2 <- textTile(x=0.8075,
                  y=0.975,
                  labels="shifts these budgets",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)
                  
vertmark <- linesTile(x=c(1,1), y=c(0,1), plot=1)

limits <- c(0.75,1.175)
at <- c(0.85,0.9,0.95,1,1.05,1.1,1.15)
labels <- c("-15%","-10%","-5%", "0%", "+5%", "+10%","+15%")

# Create plot and save to pdf -- first pass
file <- paste0("demogBudgetEffects", specname)
tc <- tile(collectedTraces,
     vertmark, text1, text2,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file, width=6.5),
     special=list(fontsizeRopeladder=5.5)
     )

ycoords <- NULL
for (i in 1:length(lhs)) ycoords <- rbind(ycoords, tc$traces[[i]]$y)

peALL0 <- peALL
pch <- peALL0*0 + 16
col <- matrix("white", nrow=nrow(peALL), ncol=ncol(peALL))
for (i in 1:ncol(peALL)) {
    for (j in 1:nrow(peALL)) {
        sig95 <- ((lwALL[j,i]-1)*(upALL[j,i]-1))>0
        sig90 <- ((lwALL0[j,i]-1)*(upALL0[j,i]-1))>0
        
        if (sig95&sig90) {
            peALL0[j,i] <- NA
        } else {
            if (sig90) {
                col[j,i] <- "gray65"
            }
        }
    }
}

nsTrace <- pointsTile(x=as.numeric(t(peALL0)),
                      y=as.numeric(ycoords),
                      col=as.character(t(col)),
                      pch=pch,
                      alpha=1,
                      layer=1,
                      size=0.51,
                      plot=1)

# Create plot and save to pdf -- second pass
file <- paste0("demogBudgetEffects", specname)
tile(collectedTraces, nsTrace,
     vertmark, text1, text2,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file, width=6.5),
     special=list(fontsizeRopeladder=5.5)
     )



## Make Plot 6

peALL <- simfrestall[14:17,]
lwALL <- simfrlowall[14:17,]
upALL <- simfrupall[14:17,]
lwALL0 <- simfrlowall0[14:17,]
upALL0 <- simfrupall0[14:17,]

colnames(peALL) <- lhsWide

row.names(peALL) <- c("Northeast",
                      "Midwest",
                      "South",
                      "West"
                      )

collectedTraces <- collectedBlanks <- vector("list", nrow(peALL))
for (i in 1:ncol(peALL)) {

    labelscur <- colnames(peALL)
    
    pch <- rep(1,nrow(peALL))
    for (j in 1:nrow(peALL))
        if (((lwALL[j,i]-1)*(upALL[j,i]-1))>0) {
            pch[j] <- 16
        }
    
    collectedTraces[[i]] <-  ropeladder(x=peALL[,i],
                                        lower=lwALL[,i],
                                        upper=upALL[,i],
                                        labels=row.names(peALL),
                                        sublabels=labelscur[i],
                                        sublabelsX=0.135,
                                        entryheight=0.375,
                                        subentryheight=1.10,
                                        size=0.6,
                                        pch=pch,
                                        spaceAbove=0.025,
                                        plot=1
                                        )
}

text1 <- textTile(x=0.778,
                  y=0.975,
                  labels="Region compared to the US as a whole",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)

text2 <- textTile(x=0.8075,
                  y=0.975,
                  labels="",
                  clip="off",
                  fontface=4,
                  cex=0.75,
                  plot=1)
                  
vertmark <- linesTile(x=c(1,1), y=c(0,1), plot=1)

limits <- c(0.75,1.175)
at <- c(0.85,0.9,0.95,1,1.05,1.1,1.15)
labels <- c("-15%","-10%","-5%", "0%", "+5%", "+10%","+15%")

# Create plot and save to pdf -- pass 1
file <- paste0("regionBudgetEffects",specname)
tc <- tile(collectedTraces,
     vertmark, text1, text2,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     yaxis = list(add=FALSE),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file, width=6.5),
     special=list(fontsizeRopeladder=5.5)
     )

ycoords <- NULL
for (i in 1:length(lhs)) ycoords <- rbind(ycoords, tc$traces[[i]]$y)

peALL0 <- peALL
pch <- peALL0*0 + 16
col <- matrix("white", nrow=nrow(peALL), ncol=ncol(peALL))
for (i in 1:ncol(peALL)) {
    for (j in 1:nrow(peALL)) {
        sig95 <- ((lwALL[j,i]-1)*(upALL[j,i]-1))>0
        sig90 <- ((lwALL0[j,i]-1)*(upALL0[j,i]-1))>0
        
        if (sig95&sig90) {
            peALL0[j,i] <- NA
        } else {
            if (sig90) {
                col[j,i] <- "gray65"
            }
        }
    }
}

nsTrace <- pointsTile(x=as.numeric(t(peALL0)),
                      y=as.numeric(ycoords),
                      col=as.character(t(col)),
                      pch=pch,
                      alpha=1,
                      layer=1,
                      size=0.51,
                      plot=1)

# Create plot and save to pdf -- pass 2
file <- paste0("regionBudgetEffects",specname)
tile(collectedTraces, nsTrace,
     vertmark, text1, text2,
     limits = limits,
     width=list(null=4),
     height=list(topaxistitle=2.5, xaxistitle=2.5),
     xaxis= list(at = at, labels = labels,cex=0.8),
     yaxis = list(add=FALSE),
     topaxis= list(at = at, labels = labels, add = TRUE,cex=0.8),
     topaxistitle = list(labels="Cumulative percent change in budget after 4 years"),
     gridlines=list(type="xt"),
     output=list(outfile=file, width=6.5),
     special=list(fontsizeRopeladder=5.5)
     )
